/* Oh, hello there. You must be looking for the file that was used to generate 
   the analysis for <DOI GOES HERE>, publised at JPART. If so, then you are in
   the right place. If not, well, better luck next time. */

/*--------Install necessary packages and schemes if not already present-------*/

// Add any required ado files within the quotations separated by a single space
local required_ados "putexcel estout" 

// Add any required schemes within the quotations separated by a single space
local required_schemes "burd"

foreach x of local required_ados {
	capture findfile `x'.ado
	// Check to see if the returned error code is equivalent to 'not found' (601)
	if _rc == 601 {
		// If it's not present, install the .ado
		ssc install `x'
	}
} 
 
foreach x of local required_schemes {
	capture set scheme `x'
	//Check to see if the returned error code is equivalent to 'not found' (111)
	if _rc == 111 {
		//If it's not present, install the scheme
		ssc install scheme-`x'
	}
} 

/*------------------------------Set global macros-----------------------------*/

// User's directory structure to get to the folder where the project data are 
// located. Change this as necessary

global path "C:\Users\mallory.compton\Dropbox\UI Payment Projects (Shared Folder)\Administrative Errors and Diversity\JPART\reptest"


/*---------------------------Clear Memory and Setup---------------------------*/

clear
set more off
label drop _all
set scheme burd

// Load analysis-ready file
use "$path\AED_JPART_replication_data.dta", clear

// Declare the data as year-interval panel data
xtset claimid claim_year, yearly


/*------------------------------Generate Figure 1-----------------------------*/

// White total
bysort claim_year: egen wh_tot = count(wkmode) if r_eth == 1 & wkmode != 999
// Black total
bysort claim_year: egen bk_tot = count(wkmode) if r_eth == 2 & wkmode != 999
// Hispanic total
bysort claim_year: egen hs_tot = count(wkmode) if r_eth == 5 & wkmode != 999
// API total
bysort claim_year: egen ap_tot = count(wkmode) if r_eth == 3 & wkmode != 999
// AIAN total
bysort claim_year: egen ai_tot = count(wkmode) if r_eth == 4 & wkmode != 999

// White by telephone	   
bysort claim_year: egen wh_tel = count(wkmode) if r_eth == 1 & wkmode == 1
// White by in-person
bysort claim_year: egen wh_per = count(wkmode) if r_eth == 1 & wkmode == 2
// White by internet
bysort claim_year: egen wh_net = count(wkmode) if r_eth == 1 & wkmode == 3
// White proportion telephone
bysort claim_year: gen wh_proptel = wh_tel/wh_tot
// White proportion in-person
bysort claim_year: gen wh_propper = wh_per/wh_tot
// White proportion internet
bysort claim_year: gen wh_propnet = wh_net/wh_tot

// Black by telephone
bysort claim_year: egen bk_tel = count(wkmode) if r_eth == 2 & wkmode == 1
// Black by in-person
bysort claim_year: egen bk_per = count(wkmode) if r_eth == 2 & wkmode == 2
// Black by internet
bysort claim_year: egen bk_net = count(wkmode) if r_eth == 2 & wkmode == 3
// Black proportion telephone
bysort claim_year: gen bk_proptel = bk_tel/bk_tot
// Black proportion in person
bysort claim_year: gen bk_propper = bk_per/bk_tot
// Black proportion internet
bysort claim_year: gen bk_propnet = bk_net/bk_tot

// Hispanic by telephone
bysort claim_year: egen hs_tel = count(wkmode) if r_eth == 5 & wkmode == 1
// Hispanic by in-person
bysort claim_year: egen hs_per = count(wkmode) if r_eth == 5 & wkmode == 2
// Hispanic by internet
bysort claim_year: egen hs_net = count(wkmode) if r_eth == 5 & wkmode == 3
// Hispanic proportion telephone
bysort claim_year: gen hs_proptel = hs_tel/hs_tot
// Hispanic proportion in person
bysort claim_year: gen hs_propper = hs_per/hs_tot
// Hispanic proportion internet
bysort claim_year: gen hs_propnet = hs_net/hs_tot

// API by telephone
bysort claim_year: egen ap_tel = count(wkmode) if r_eth == 3 & wkmode == 1
// API by in-person
bysort claim_year: egen ap_per = count(wkmode) if r_eth == 3 & wkmode == 2
// API by internet
bysort claim_year: egen ap_net = count(wkmode) if r_eth == 3 & wkmode == 3
// API proportion telephone
bysort claim_year: gen ap_proptel = ap_tel/ap_tot
// API proportion in person
bysort claim_year: gen ap_propper = ap_per/ap_tot
// API proportion internet
bysort claim_year: gen ap_propnet = ap_net/ap_tot

// AIAN by telephone
bysort claim_year: egen ai_tel = count(wkmode) if r_eth == 4 & wkmode == 1
// AIAN by in-person
bysort claim_year: egen ai_per = count(wkmode) if r_eth == 4 & wkmode == 2
// AIAN by internet
bysort claim_year: egen ai_net = count(wkmode) if r_eth == 4 & wkmode == 3
// AIAN proportion telephone
bysort claim_year: gen ai_proptel = ai_tel/ai_tot
// AIAN proportion in person
bysort claim_year: gen ai_propper = ai_per/ai_tot
// AIAN proportion internet
bysort claim_year: gen ai_propnet = ai_net/ai_tot

// Hoo boy that was a lot of work just to make one line graph eh? Oh well!   	  
line wh_proptel wh_propper wh_propnet bk_proptel bk_propper bk_propnet hs_proptel hs_propper hs_propnet claim_year, ///
       legend( ///
	           size(small) ///
			   label( 1 "White: Telephone" ) ///
			   label( 2 "White: In-Person" ) ///
			   label( 3 "White: Web" ) ///
			   label( 4 "Black: Telephone" ) ///
			   label( 5 "Black: In-Person" ) ///
			   label( 6 "Black: Web" ) ///
			   label( 7 "Hispanic: Telephone" ) ///
			   label( 8 "Hispanic: In-Person" ) ///
			   label( 9 "Hispanic: Web" ) ///
			   col(3) ///
			   position(12) ///
			 ) ///
	   xlabel(2000(2)2018) ///
	   xtitle("Claim Year") ///
	   ytitle("Proportion of Total Claims Filed")
graph export "$path\Figure_1_Method_of_claim_submission_by_race_ethnicity_over_time.tif", as(tif) replace


/*------------------------------Generate Analysis-----------------------------*/

// Race/Ethnicity: r_eth : White Non-Hispanic as base case
// Filing method: icfilmeth OR wkfilmeth : telephone as base case
// Method auditor contacted claimant: methinfoobt : Mail or other method as base case
// Gender: female : Male as base case
// Age: cage : continuous (try base and polynomial separately)
// Education: edu : HS as base case
// Training: voctech : Never attended as base case
// Nationality: citizen : US citizen as base case
// Reservation wage: lnreswg : continuous (logged)
// Prior wages: lnbpw : continuous (logged)
// Union status: union: Not union as base case
// Whether dependents were claimed: deps : No dependends claimed as base case
// Program Type: altprog : Standard UI as base case
// Work search requirements: wsreq : Required as base case
// Labor Exchange registration required: labex : Requred as base case
// State Fixed Effects: st_fips : CA as base case
// Year Fixed Effects: 2005 as base case
// Quarter Fixed Effects: Q2 as base case
// Employer Industry Code: Manufacturing as base case

// Control variable macros
global indiv_controls "i.female c.cage c.agesq ib2.edu i.citizen"
global policy_controls "i.altprog"
global state_fe "st_fips"
global year_fe "claim_year"
global quart_fe "claim_quart"
global indus_fe "empcode"

//Clear saved estimates
eststo clear


/*-----------------------------Estimate Model 1-------------------------------*/

log using "$path\AED_JPART_Model_1.smcl", replace	

eststo m1: logit agfault ///
      i.r_eth##i.wkmode ib3.methinfoobt ///
	  $indiv_controls ///
	  $policy_controls ///
	  ib06.$state_fe ib2005.$year_fe ib2.$quart_fe ib31.$indus_fe ///
	  if wkmode != 999 ///
	  , vce(cluster st_fips) or
   
log close

// Generate Figure 2  
margins r_eth, over(wkmode)
marginsplot, ///
	recast(scatter) ///
	level(83) ///	
	ytitle("Predicted Probability of Error", size(medium)) ///
	yscale(range(0(0.025)0.2) extend) ///
	legend(ring(0) position(11) bmargin(0 20 20 0)) ///
	title("Any Agency Error", size(medium) position(12)) ///
	name(mod1, replace)	  
graph export "$path\Figure_2_AME_any_agency_error_by_race_ethnicity.tif", as(tif) replace
  

/*----------------------------Generate Descriptives---------------------------*/

// NOTE: Run descriptives after Model 1 to condition on inclusion in sample

log using "$path\AED_JPART_descriptives.smcl", replace
	
foreach var of varlist agfault r_eth wkmode methinfoobt female edu citizen altprog {
	tab `var' if !missing(agfault, cage, r_eth, wkmode, methinfoobt, female, /// 
	                      edu, citizen, altprog) ///
	             & wkmode != 999 ///
				 & e(sample)
}

sum cage if !missing(agfault, cage, r_eth, wkmode, methinfoobt, female, edu, ///
                     citizen, altprog) ///
	             & wkmode != 999 ///
				 & e(sample)	

log close	  
	  
	  
/*-----------------------------Estimate Model 2-------------------------------*/

log using "$path\AED_JPART_Model_2.smcl", replace	

eststo m2: mlogit c_agentnomial ///
      i.r_eth##i.wkmode ib3.methinfoobt ///
	  $indiv_controls ///
	  $policy_controls ///
	  ib06.$state_fe ib2005.$year_fe ib2.$quart_fe ib31.$indus_fe /// 
	  if wkmode != 999 ///
	  & hsamp == 1 ///
	  , vce(cluster st_fips) rrr

log close


// Generate Figure 3. I recommend stepping away from the computer to perhaps
// enjoy a delicious beverage. Or maybe do some yoga. Either way, unless you're
// running this on a high-performance cluster system (and if you are, kudos to
// you for having a University that lets you do that kind of thing), this is
// going to take a long time.

	  
// MFX: Race and Ethnicity by Claim Filing Technology
margins r_eth, over(wkmode) predict(outcome(1))
marginsplot, ///
	recast(scatter) ///
	level(83) ///
	ytitle("Predicted Probability of Error", size(medium)) ///
	legend(ring(0) position(11) bmargin(0 20 20 0)) ///
	title("Overpayment", size(medium) position(12)) ///
	name(m2_out_1, replace)
//graph export "$path\m2_mfx_out1.png", as(png) replace	

margins r_eth, over(wkmode) predict(outcome(2))
marginsplot, ///
	recast(scatter) ///
	level(83) ///
	ytitle("") ///
	legend(ring(0) position(11) bmargin(0 20 20 0)) ///
	title("Underpayment", size(medium) position(12)) ///
	name(m2_out_2, replace)
//graph export "$path\m2_mfx_out2.png", as(png) replace

margins r_eth, over(wkmode) predict(outcome(3))
marginsplot, ///
	recast(scatter) ///
	level(83) ///
	ytitle("") ///
	legend(ring(0) position(11) bmargin(0 20 20 0)) ///
	title("Denial, All Other", size(medium) position(12)) ///
	name(m2_out_3, replace)
//graph export "$path\m2_mfx_out3.png", as(png) replace	

margins r_eth, over(wkmode) predict(outcome(4))
marginsplot, ///
	recast(scatter) ///
	level(83) ///
	ytitle("Predicted Probability of Error", size(medium)) ///
	legend(ring(0) position(11) bmargin(0 20 20 0)) ///
	title("Denial, Non-Mon/Sep", size(medium) position(12)) ///
	name(m2_out_4, replace)
//graph export "$path\m2_mfx_out4.png", as(png) replace

* Combining primary outcomes into single small multiple graph
graph combine m2_out_1 m2_out_2 m2_out_4 m2_out_3 , col(2) ycommon
graph export "$path\Figure_3_AME_Agency_Errors_By_Type_By_Race_Ethnicity_By_Technology.tif", ///
             as(tif) replace width(1950) height(1650)
	  
  
/*------------------------OUTPUT ANALAYSIS TO TABLES--------------------------*/
	  
	  
// Output the models to csv
esttab m1 using "$path\model_1_results.csv", ///
	   cells("b(star fmt(3)) se ci_l ci_u") ///
	   eform ///
	   stats(N r2_p) /// 
	   starlevels(* 0.05 ** 0.01 *** 0.001) ///
	   stard ///
	   indicate("State Fixed Effects = *.$state_fe" ///
	            "Year Fixed Effects = *.$year_fe" ///
				"Quarter Fixed Effects = *.$quart_fe" ///
				"Industry Fixed Effects = *.$indus_fe") ///
	   wide plain label ///
	   replace 
	   
esttab m2 using "$path\model_2_results.csv", ///
	   cells("b(star fmt(3)) se ci_l ci_u") ///
	   eform ///
	   stats(N r2_p) /// 
	   starlevels(* 0.05 ** 0.01 *** 0.001) ///
	   stard ///
	   indicate("State Fixed Effects = *.$state_fe" ///
	            "Year Fixed Effects = *.$year_fe" ///
				"Quarter Fixed Effects = *.$quart_fe" ///
				"Industry Fixed Effects = *.$indus_fe") ///
	   wide plain label ///
	   replace 


/*----------------------------THAT'S IT THAT'S ALL----------------------------*/
